% Solve for baseline U.S. economy

baseline.tau    = pr.tau;
baseline.lambda = fzero(@(x) BaselineAllocation(x,pr,w,pai),1);
BaselineAllocation(baseline.lambda,pr,w,pai);

load c
load y
load G
baseline.c       = c;
baseline.y       = y;
baseline.output  = sum(pai.*y);
baseline.welfare = sum(pw.*pai.*Utility(c,y,w,pr));
baseline.MTR     = 1-pr.Om./c.^(-pr.gam).*w.^-(1+pr.sig).*y.^pr.sig;
baseline.ATR     = (y-c)./y;
baseline.EfAveMTR= sum(pai.*baseline.MTR.*baseline.y/baseline.output);
[temp.y,temp.i]  = min(baseline.y);
temp.c           = baseline.c(temp.i);
baseline.Tr2Y    = (temp.c-temp.y)./baseline.output;
baseline.TrG2Y   = (temp.c-temp.y+G)./baseline.output;

save ./variables/baseline baseline
delete c.mat y.mat 


% Pareto on family consumption and income
if (index_case>=3)
    coef_c = regress(log(pai(end-100:end)),[0*w(end-100:end)+1 log(baseline.c(end-100:end))]);
    coef_y = regress(log(pai(end-100:end)),[0*w(end-100:end)+1 log(baseline.y(end-100:end))]);
    %display(['Pareto on income, consumption = ',num2str(-coef_y(2)),', ',num2str(-coef_c(2))]) % ratio is 1/(1-tau)
    lambda_c = -coef_c(2);
else
    lambda_c = Inf;
end


cumsum_pai = cumsum(pai);
deciles = ones(10,1);
a = 1; 
for i = 1:9
    b = find(cumsum_pai>0.1*i,1,'first');
    deciles(i) = sum(pai(a:b).*baseline.y(a:b))/sum(pai.*baseline.y);
    a = b+1;
end

deciles(10) = sum(pai(a:end).*baseline.y(a:end))/sum(pai.*baseline.y);
quintiles = ones(5,1);
for i  =1:5
    quintiles(i,1) = sum(deciles(2*(i-1)+1:2*i));
end
